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1.  Introduction. 


The  problem  of  computing  the  variance  of*  ample  of  N  data  potato  {«,)  ie  < 
at  first  glance,  to  be  slmoet  trivial  but  can  in  tact  be  quite  dURcult,  particularly 
and  tbe  variance  Is  small.  The  fundamental  ratrilallen  caorisU  ef  computing  tb 
of  the  deviations  from  tbe  mean, 
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The  sample  variance  is  then  SfN  or  8/(N  - 1)  depending  on  the  application.  The  formulas  (1.1) 
define  a  straightforward  algorithm  for  computing  8.  This  will  be  called  the  standard  two  pam 
algorithm,  since  It  requires  passing  through  the  data  twice:  once  to  compute  t  and  then  sqain  to 
compute  5.  This  may  be  undesirable  in  many  applications,  for  example  when  the  data  sample  Is 
too  large  to  be  stored  in  main  memory  or  when  the  variance  Is  to  be  dynamically  as  the 

data  is  collected. 

In  order  to  avoid  the  two  pees  nature  of  (1.1),  H  Is  standard  practice  to  the 

definition  of  S  into  the  form 

*-£•?-*(£«.)*  («•») 


This  form  is  frequently  suggested  in  statistical  textbooks  and  will  be  —n— <  the  rat  jraw 

algorithm.  Unfortunately,  although  (14)  b  mathematically  equivalent  to  (1.1),  numerically  it  can 
be  disastrous.  The  quantities  £*?and  $(£*)•  may  be  very  large  in  practice,  and  will  generally 
be  computed  with  some  rounding  error.  If  the  variance  b  small,  these  numbers  should  >i  out 
almost  completely  in  the  subtraction  of  (14).  Many  (or  aD)  or  the  correctly  compute  digits  will 
cancel,  leaving  a  computed  5  with  a  possibly  unacceptable  relative  error.  Tbe  computed  S  can 
even  be  negative,  a  blessing  in  disguise  since  thb  at  bast  alerts  the  programmer  that  disastrous 
cancellation  has  occured. 

To  avoid  these  difficulties,  several  alternative  one-pass  algorithms  have  been  Introduced.  These 
Include  the  updating  algorithma  of  Youngs  and  Cramer(Il),  Wdford[IO),  Wcot|fi],  Uanson[t],  and 
CottonJfiJ,  and  the  pairwise  algorithm  of  the  present  authors!!].  In  describing  these  algorithms  we 
will  use  the  notation  T,,  and  My  to  denote  the  sum  and  the  mean  of  the  data  points  a,  through 
*i  respectively, 


**  “  (j-i+l)7* 


and  Sij  to  denote  the  sum  of  squares 


*#  -  £<«»  - 


«k)*. 


For  competing  an  unweighted  sum  of  squares,  as  we  consider  here,  the  algorithms  oT  Watford,  West 
and  Hanson  are  virtually  identical  and  are  based  on  the  updating  formulas 

hitj  *  hdij-i  ♦  j(*y  -  hiij-i) 


*u  -  8,J-,  *u- «X*< -  Hi 


1 


(I4a) 

(14*) 


with  Mi.i  *  si  and  S|,i  *  0.  The  ddwd  value  of  8  b  ultimately  obtained  aa  8%^.  Tbs 
updating  formulas  of  Youop  and  Cramer  are  dmtltr. 

TU  “  TU-t  ♦  */  M 

s,  J  -  Su.t  +  -  ru)*  («-•»> 

with  Tt,\  ■»  si  and  S\,\  “  0.  These  two  algorithms  have  dollar  numerical  behavior  aad  are 
more  stable  than  the  textbook  algorithm.  Note,  la  particular,  that  with  both  of  these  algorithms 
$  mt  Si'N  la  computed  aa  the  sum  of  aonnegative  quantities.  Cotton's  update  la  no  more  stable 
than  the  textbook  algorithm  and  should  not  be  aaod  (see  [1]). 

The  updating  formulas  (1.4)  can  be  generalised  to  allow  as  to  eomblne  two  samples  of  arbitrary 
else.  Suppose  we  have  two  samples  {*}£.!,  {sf)£22+i  and  we  know 

*<»  J  q* 

4m  1  4mm Ml 

8l,m  *  ]^(*»  —  "“I’ll*)* i  “  J)  (•<  "  * 

iml  "  fc— <1 

Then,  If  we  combine  all  of  the  data  Into  a  anmple  of  rise  m  ♦  a,  we  have 

Ti,m+m  *  1U  +  Vfl«a4« 

8i,m+n  •  d||H  4* 

-)’•  <14*> 

When  m«s  this  reduces  to 

8l,tm  *  8i,m  +  +  jJJjllVe  "  la+lifc*)* •  0*®) 

This  formula  forms  the  bads  of  the  pairwise  algorithm.  The  pairwise  summation  algorithm  for 
computing  the  sum  of  N  numbers  is  wall  known  aad  can  be  described  recursively  by  stating  that 
Ti.tm  shall  be  computed  as 

?ttSm  ■  Ji,u  + 

with  each  of  the  sums  on  the  right  band  ride  computed  la  a  rimilar  manner.  Formula  (1.6) 
defines  the  analogous  pairwise  algorithm  for  computing  the  variance.  This  can  be  implemented 
In  a  one* pass  manner  using  only  0(log  Af)  internal  storage  locations  as  discussed  in  (2]  and  also 
by  Nash[7).  All  logarithms  in  this  paper  are  base  S.  It  can  easily  be  shown  that  the  use  of  the 
pairwise  summation  algorithm  reduces  relative  errors  in  T|^  from  O(N)  to  0(\ogN)  tm  N  -*ao. 
The  pairwise  variance  algorithm  can  be  mpeeted  to  have  the  came  advantage,  as  Is  confirmed 
numerically. 

Incidentally,  palrwice  summation  cna  be  used  in  Implementing  (1.1)  (both  in  computing  »  aad 
In  forming  8)  or  (1 J)  with  rimilar  benefits. 

Other  devices  can  else  be  need  to  increase  the  accureey  of  the  computed  S.  Fsr  data  with  a 
large  mean  value  t,  experience  has  shewn  that  substantial  galas  in  accuracy  can  be  achieved  by 
shifting  all  of  the  data  by  same  approximation  to  •  before  attempting  to  compute  5.  Even  n  crude 
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estimate  of  t  can  yield  dramatic  Hnpror— eats  hi  accuracy,  m  wu  Mod  Ml  smart  to  •  tass-paas 
algorithm  la  order  to  first  estimate  8.  This  is  dbcassod  la  detail  b  metiee  S.  Dosws,  whea  the 
shift  e»  the  eompatod  mean  aad  the  tart  hoot  algorithm  (14)  b  thea  applad  to  tits  drifted  data, 
om  obtains  the  oorracfod  two-pom  dprttta 

■  CW) 

Hera  the  first  term  is  rimply  the  Im  psm  algorithm  (1.1a).  The  msaad  tana  weald  ho  amo  la 
exact  computation,  hut  In  practice  is  a  good  apprmrtmsHoa  to  the  error  la  the  first  tana.  Nolo 
that  In  this  case  use  of  the  textbook  algorithm  dam  aot  load  to  cataatrophk  caacsUallon,  riaee 
the  correction  Is  generally  much  smaller  than  the  first  term.  This  algorithm  was  first  printed  oat 
to  the  authors  by  Professor  A.  B98rck(l]  who  amigaalail  this  enrradlea  term  booed  solely  aa  the 
error  aaalyris  of  the  two- pass  algorithm^].  4a  alternative  (aad  impruved)  error  aaalyris  Is  riven 
la  section  9. 

Initially  algorithms  for  computing  the  eerie  nee  ware  judged  solely  aa  the  haste  ef  empirical 
studies(f],  [•],  (11).  More  recently  rigorous  error  heaade  hews  hooa  obtained  for  nmay  algorithms)!). 
(I),  (4).  Our  aim  here  is  to  present  a  aaifiad  eanrey  of  error  aaalysao  for  the  eh  era  awatinood  al¬ 
gorithms  and  techniques.  Some  of  this  material  ie  hefieead  to  ha  now,  particularly  the  laeaaHgetina 
Into  the  effects  of  shifting  the  data.  Based  aa  this  sarvuy,  specific  reonmmaadeWnas  will  ho  amde 
as  to  which  algorithm  should  be  aeed  la  aariaas  aoataxta. 

S.  Condition  numbers  and  error  naalyrio 

Chan  and  Lswis)4)  first  derivad  the  eamdltiea  number,  «,  of  a  eemple  {*J  (with  raapoct  to 
computing  the  variance).  This  condition  number  meaeuree  the  aandtirity  of  S  for  the  given  data 
set.  If  relative  errors  of  else  7  are  introduced  into  the  a*,  than  the  relative  change  in  S  is  bounded 
by  07.  Chan  and  Lewis  showed  this  to  be  true  ap  to  0(7*).  In  fact  it  is  strictly  true  as  noted  by 
van  Nes(8).  Physical  data  almost  always  has  some  uncertainty  fee  It,  aad  this  uncertainty  will  be 
magnified  by  the  factor  a  in  5.  If  nothing  rise,  errors  are  introduced  In  representing  the  data  on 
the  computer,  mod  so  n  value  of  S  computed  on  n  computer  with  machine  accuracy  u  may  have 
relative  errors  as  large  aa  mi  ragardlem  of  what  algorithm  ie  oeed.  This  value  mi  can  be  need  as  a 
yardstick  by  which  to  judge  the  accuracy  of  tho  various  algorithms,  mpodolly  Since  error  bounds 
can  often  be  derived  which  are  functions  solely  of  s,  u,  aad  N, 

If  we  define  the  2- norm  ef  the  data  by 


M8-Z.J. 

4*1 

then  the  condition  number  for  this  problem  is  given  by 

..if. 

vs 

When  8  io  small  aad  t  is  not  dose  to  aero  we  obtain  the  useful  approximation 

« m  9^/nJS  (for  S  wnall,  t  aonsero)  (14) 

which  is  tho  moan  divided  by  tho  standard  deviation.  We  always  hove  *  £  1,  and  in  many 
ritueUons  a  is  very  large. 


Vn/S. 


(*•») 


Table  2.1  shows  the  asymptotic  error  Woods  for  the  algorithms  discerned.  These  ore  Woods 
on  the  relative  error  |(S  -  S)/S\  io  the  eompwted  voles  S.  Small  ssostsot  msltipDsrs  Wes  Wsa 
dropped,  for  clarity.  Higher  order  terms  have  oho  Wsa  dropped,  Wt  the  terms  shown  dsmieete 
the  error  Woods  whenever  the  relative  error  h  leas  thao  1.  The  Woods  for  the  textbook  algorithm 
and  West’s  updating  are  derived  by  Chan  aad  Lowta(4].  TW  two  pees  error  Wood  lododlog  tW 
term  (which  can  dominate  lo  practice)  is  derived  la  (3).  Bo  ends  for  these  algorithms  oaiag 
pairwise  summation  can  be  found  similarly.  The  pairwise  variance  algorithm  booed  is  e  conjecture 
Weed  on  the  form  of  the  error  Wund  for  Youngs  and  Cramer  epdatiag  and  nrperimeotal  results. 
Tbe  error  analysis  for  the  corrected  two-paa  algorithm  is  given  la  aectioa  3. 

Graphs  of  these  Wunds  are  shown  h  Figures  3.1.  through  24  along  with  some  taper imeotai 
results.  Each  plot  has  a  on  the  abockaa  and  the  relative  error  ia  8  oo  the  ordinate.  TW  lower  carve 
in  each  figure  shows  the  error  Wund  for  N  ■■  t4,  the  upper  com  4M6.  TW  oumorieal 

experiments  were  performed  on  an  IBM  S081  computer  at  the  Stanford  linear  Accelerator  Center. 
The  data  used  was  provided  by  a  normal  random  number  generator  with  mean  1  aad  a  variety  of 
different  variances  1  ^  e*  £  10"1*.  Fbr  this  choice  of  the  mens,  a  m  1/e  (see  (24)).  In  onch 
case  the  results  have  been  averaged  over  SO  rone.  Slagle  precision  was  need  ia  all  of  tW  taste, 
with  machine  accuracy  mix  10-T.  TW  “eerreet"  onswer  fbr  one  in  computing  the  error  woe 
calculated  in  double  precision.  Tbs  moulting  arrow  are  denoted  In  tbe  figures  by  the  symbols  + 
(for  N  mm  04)  and  X  (for  AT  -  40M). 

Tbe  experimental  results  confirm  tbe  general  form  of  the  error  Woods  given  In  Thble  2.1. 
In  particular  the  graph*  for  tbe  two- peas  algorithms  shew  Ww  the  higher  order  terms  (such  aa 
N* **u*)  begin  to  dominate  the  error  at  fairly  modest  values  of  a 


Table  3.1.  Error  Wunds  for  tbe  relative  error  |*§®|  ^  *k«  computed  value  S.  Only  tbe  dominant 
terms  are  shown,  and  small  constant  factors  have  Wen  suppressed  for  clarity. 


1.  textbook 

ff**s 

3.  textbook  with  pairwise  summation 

n*ulogN 

3.  two- pass 

N*  +  AT**  V 

4.  two-pass  with  pairwise  summation 

o  log  N  +  (c*  log  N)* 

5.  corrected  two- pass 

No  +  N***** 

0.  corrected  two-pass  with 

pairwise  summation 

n  log  N  +  c*«*  log*  N 

7.  updating 

N* « 

S.  pairwise 

no  log  N  (conjectured) 
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3.  Computation*  with  shifted  data 

IT  we  replace  the  original  data  {a,}  by  dlfted  data  {!»}  Mud  by 

(U) 

for  some  fixod  shift  d,  then  the  new  data  baa  mu  t  —  4  aad  8  remains  unchanged  (aaaamiag 
tbe  it  are  computed  exactly).  la  practice  data  with  a  nans  ere  mean  la  frequently  ablftod  by  acme 
a  priori  estimate  of  tbe  mean  before  attempting  to  aampete  5.  Thb  will  generally  Increase  tbe 
accuracy  of  tbe  computed  S.  We  will  analyse  this  improvement  by  investigating  tbe  depen  dee  re 
of  tbe  condition  number  on  tbe  shift.  Bound#  on  £, tbe  condition  nambtr  of  tbe  ekiftod  data,  arc 
derived  for  various  choices  of  tbe  shift  4.  Tbsse  can  than  be  inserted  in  place  of  « la  tbe  bounds 
of  Tabic  3.1  to  obtain  error  bounds  for  each  of  tbe  algorithms  with  drifted  data. 

From  the  definition  of  tbe  condition  number  we  have 

*•-«+ (mi 

Comparing  this  with  (2.1)  we  eee  that  £  <  «  whenever  |d  - 1|  <  ft|,  La.,  whenever  4  Bee  between 
0  and  22.  Taking  rfsl  gives  perfectly  conditioned  data,  i*l.  h  practice  we  caaaot  compute  g 
exactly  and  usually  will  not  even  attempt  to  compute  H  (except  when  using  a  tro-psa  algorithm). 
Instead,  we  use  some  rough  estimate  which  is  easily  computed  without  a  separate  pass  through  all 
of  the  data. 

Frequently  a  shift  4  is  obtained  by  timply  "eyeballing”  tbe  data.  8ueh  a  technique  aright  be 
expected  to  yield  an  approximation  4  which  is  within  a  few  standard  deviations  of  the  mean.  This 
is  sufficient  to  give  completely  satisfactory  bounds  on  £.  Recall  that  the  standard  deviation  b 
y/S/N  and  suppose  that  |f  -  4[  <  p y/$/N  for  asms  small  p.  Then  (3.2)  gives 

l^l+F*.  (U) 

For  example,  if  d  is  within  one  standard  deviation  of  tbe  mean  then  £  <  >/$■  This  result  b 
completely  independent  of  5  aad  N. 

It  b  not  always  possible  to  obtain  an  approximation  in  thb  manner,  nor  b  it  always  valid  to 
make  such  an  assumption  on  its  accuracy.  Another  bound  on  £  can  be  easily  obtained  by  aasuming 
only  that 

■in**  £  *£ ■»**.•• 

This  b  easily  guaranteed,  for  etample  by  choosing  one  of  the  data  points  as  the  shift.  When 
min  *i<4<  max  **,  we  have  {*  -  d)*  <  £^(X  -  *<)*  ■*  8  and  so  from  (3.3), 

£*  <  1  +  AT.  (3.4) 

Thb  bound  b  not  as  mtbfactory  as  (3.3),  but  for  moderate  values  of  N  it  may  be  sufficient  to 
guarantee  acceptable  errors  in  S. 

For  tbe  esse  in  which  we  shift  by  s  tingle  data  point,  4  -«  x,  for  aome  j,  we  can  obtain  some 
interesting  probabilistic  refinements  of  (3.4).  Equality  in  (3.4)  b  unattainable  and  approximate 
equality  bolds  only  when  ^  ^ 

i 

Ls.,  only  when  i,  Bes  considerably  farther  from  X  than  do  any  of  the  other  x*.  If  x,  b  picked  at 
random  from  tbe  sample  {«<},  then  the  expected  value  of  £*  will  be  much  smaller  than  1  +  N.  In 
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fact,  aoce  E|(*  -  *.)*] «  S/S,  (tb.  MbIUm  of  (bo  aapb  wtaaoo),  m  bon  from  (U)  Ibok 

*P*l-«  («) 

Independent  of  AT  and  S.  Note  that  this  In  alno  Independent  of  the  underlying  dletiibetion  at  tke 
{*%}.  We  aenomed  only  that  *,  van  cheoen  from  {«.}  vtth  a  uniform  dintribattoa.  Aftsrnattvety 
ve  could  ehoone  the  daU  value  with  a  faced  Index,  eny  and  ansa  me  that  the  data  In  ordered 
randomly.  Thin  may  not  be  a  valid  assumption  If,  far  Initial  transients  are  present  In  the 

data. 

Improved  upper  hounds  of  the  farm  (3.4)  eaa  alno  he  obtained  probablBatkally  which  held 
with  probability  done  v.  i  For  faced  h,  1  k  jg  Af,  the  Inequality 

can  bold  for  at  most  N/k  values  s»f  ».  Otherwise  we  would  have  £(Z  -  *,)*  >  %(kS/N)  mm  8. 
Tbu*  if  x,  ia  chosen  at  random,  there  in  a  probability  of  at  least  {N  -  N/k)/N  —  1-1  /*  that 
(S  —  Xj)*  <  kS/N.  It  follows  that 

**  <1  +  *  with  probability  at  leant  1  - 1/*  for  1  <  *  £  Af.  (3.6) 

If  TV  >  100  we  have,  for  example. 


«*  <  101  vtth  probabifity  0.00. 


This  is  agsin  independent  of  Af  and  S  when  the  shift  r,  in  ehonen  at  random  tram  the  sample. 

We  can  generalise  this  choice  of  d  by  wring  the  average  of  some  p  data  points,  p  «  Af.  This 
average  will  be  denoted  by  i,  «=  £  *,  /p,  the  sum  being  over  the  chosen  p  data  points.  We  assume 
that  p  is  sufficiently  small  that  rounding  errors  In  computing  Z,  can  be  ignored.  Specifically  this 
requires  apu  <  1.  The  condition  number  corresponding  to  this  shift  Is  bounded  by  using  Cauchy's 
inequality, 

S* -(+£(,-»,)• 


* 


Si  +  *. 

f 


For  p  mt  1  this  reduces  to  (3.4). 

We  now  consider  the  case  in  which  the  computed  mean  is  used  as  the  shift.  In  general  we 
cannot  Ignore  rounding  errors  In  computing  Z.  Instead  we  compute  some  approximate  floating 
point  value  fl(i),  given  by 

+  (M) 


c-i 


where  the  6  are  bounded  by 


161  S  Nn 


(*.») 


wbcD  the  grail  (forward)  summsUoa 
can  be  replaced  by  log//.  New  we  ■ 


s  i + j^wctfiii 

si+^ldlL- 


a»») 


Here  we  have  uaed  (2.1)  and  the  teneral  Inequality  Q(Q|  <  AfUflJ*.  where  g£||a,  «  max*  |{,|.  IMag 
(3.0)  we  can  rewrite  (3.10)  ae 

S*£  1  +  N*«V.  (3.11) 

Note  that  due  to  the  dependence  oa  it,  t fee  bound  (3.11)  may  be  wetee  than  the  bounds  obtained 
for  more  primitive  estimates  of  d.  This  reflects  rituailons  which  can  actually  occur  la  practice. 
One  can  easily  construct  examples  where  the  computed  mean  does  not  erven  lie  between  min  a*  and 
max  *,  and  hence(S  —  fl(*))*  is  larger  than  max,-(f  —  ztf.  In  this  case  one  is  better  off  shifting  by 
any  single  data  point  than  by  the  computed  mean. 

Of  course  shifting  by  the  computed  mean  may  also  he  an  undesirable  choice  from  the  standpoint 
of  efficiency,  since  it  requires  a  separate  pass  through  the  data  to  compute  fl(S).  Nonetheless,  when 
a  two-pass  algorithm  is  acceptable  and  JV*«*u9  is  small  (<  1,  say),  this  shift  followed  by  a  one-pass 
algorithm  provides  a  very  dependable  method  for  computing  S.  The  corrected  two-pass  algorithm 
(1.7)  is  of  this  form,  it  consists  of  the  textbook  algorithm  on  data  shifted  by  ■(*).  Its  error  hound 
Afu(l  +  Ar**V)  is  eerily  derived  from  (3-11)  and  the  textbook  algorithm  bounds  of  Table  2.1. 

Other  one-pass  algorithms  could  also  be  used  in  conjunction  with  a  shift  by  the  computed 
mean.  However,  if  a  good  shift  has  been  chosen  so  that  i  iw  1,  alt  one-pass  algorithms  are  essentially 
equivalent  with  a  bound  N%  (or  u  log  AT  for  algorithms  wing  pairwise  summations).  Since  the 
textbook  algorithm  is  the  moot  efficient  one- peas  algorithm  (requiring  only  N  multiplications  and 
S N  additions  as  opposed  to  IN  multiplications  and  IN  additions  for  the  updating  algorithms,  for 
example),  it  is  the  method  of  choiee  except  in  rare  instances 


The  results  of  the  previous  sections  provide  a  basis  for  making  an  intelligent  choice  of  algorithm 
for  accurately  computing  the  sample  variance.  First  we  note  that  if  a  parallel  processor  is  available, 
the  data  can  be  split  up  into  mailer  samples  and  the  sum  of  squares  computed  for  each  sample 
individually.  These  can  then  be  combined,  and  the  global  sum  or  squares  computed,  by  using  the 
updating  formulas  (1.3).  In  that  ease  the  considerations  below  apply  for  each  processor. 

There  is  one  rituation  in  which  the  textbook  algorithm  (1.3)  can  be  recommended  as  It  stands. 
If  the  data  constats  only  of  integers,  small  enough  that  no  overflows  occur,  then  (1.2)  should  he 
used  with  the  sums  computed  in  integer  arithmetic.  In  this  case  no  roundoff  errors  occur  until  the 
Anal  step  of  combining  the  two  sums,  in  which  a  division  by  N  occurs. 

For  aon-intagral  data  we  must  flrot  dselde  whether  to  use  a  ono-pam  or  a  two  psm  algorithm. 
If  all  of  the  data  lie  la  hlgh-spood  memory  and  we  are  not  iotermtod  in  dynamically  updating 
the  variance  as  now  data  is  eeUoclsd,  then  a  Iwo-paas  algorithm  ta  probably  acceptable  and  the 


II 


corrected  two- put  algorithm  (1.7)  la  recoma— dod.  It  Nh  Urge  aad  Ugh  —racy  U  seeded,  H 
may  be  worthwhile  to  use  pairwise  summation  la  hnple—tlag  this  algorithm. 

If  a  one-pass  algorithm  is  to  be  esed,  the  Int  step  Is  to  drift  the  data  as  well  as  possible, 
perhaps  by  some  2,  as  discussed  in  8eetlon  3.  New  aa  appropriate  eae  pass  algorithm  amst  he 
chosen.  We  should  first  estimate  *,  the  eoadHUa  aamber  ef  the  shifted  data,  perhaps  by  aae  ef  the 
bounds  of  8ectlon  3.  If  Nk**,  the  error  bound  for  the  textbook  algorithm,  is  at  least  as  small  aa 
the  desired  relative  accuracy,  then  the  textbook  algorithm  can  bo  aeed  oa  the  drifted  data.  If  thin 
bound  is  too  Urge,  we  should  resort  to  a  less  sOdeat  algorithm  for  mfety.  Tbs  dependence  oa  N 
can  be  reduced  by  the  use  of  pairwise  summation.  The  dependence  —  £  can  be  reduced  by  adng 
an  updating  algorithm.  The  use  of  the  pairwise  algorithm  should  reduce  both  of  them  factors. 
When  /V  is  a  power  of  t  the  pairwise  algorithm  Is  fairly  easy  to  Implement  and  requires  only  IN 
multiplications  and  4 N  additions,  which  is  better  than  the  updating  algorithms.  For  general  N 
digbtiy  more  work  (particularly  human  work)  is  required,  making  it  lens  attractive. 

The  decision  procedure  just  described  is  shown  graphically  in  Figure  S.l. 
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